SUBATECH-95-05 



MARCH 1995 



HD-TVP-94-24 



Microcanonical Treatment of Hadronizing the Quark-Gluon Plasma 

Klaus WERNER 
Institut fur Theoretische Physik, Universitdt Heidelberg, Germany 

J6rg AICHELIN * 

Laboratoire de Physique Subatomique et des Technologies Associees, 
Universite de Nantes - EMN - IN2P3/CNRS, Nantes, France 



We recently introduced a completely new way to study ultrarelativistic nuclear scattering by providing a link 
between the string model approach and a statistical description. A key issue is the microcanonical treatment 
of hadronizing individual quark matter droplets. In this paper we describe in detail the hadronization of these 
droplets according to n-body phase space, by using methods of statistical physics, i.e. constructing Markov chains 
of hadron configurations. 
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I. INTRODUCTION 

Studying nuclear collisions at ultrarelativistic energies 
(-Ecms/micleon ^> 1 GeV) is motivated mainly by the ex- 
pectation that a thermalized system of quarks and gluons 
(quark-gluon plasma) is created [1]. There are essentially 
two directions for modelling such interactions: dynamical 
and thermal approaches. The former ones refer to string 
models [2-7] or related methods [8], supplemented by 
semihard interactions at very high energies [9-12]. Here, 
a well established treatment of hadron-hadron scattering, 
based on Pomerons and AGK rules [13], is extended to 
nuclear interactions. Thermal methods [14-19] amount 
to assuming thermalization after some initial time To, 
with evolution and hadronization being mostly based on 
ideal gas assumptions. 

We recently introduced a completely new approach 
[20,21], more realistic than the string model and more 
realistic than thermal approaches, providing a link 
between the two. Based on the string model, we first de- 
termine connected regions of high energy density. These 
regions are referred to as quark matter (QM) droplets. 
Presently, a purely longitudinal expansion of the QM 
droplets is assumed. Once the energy density falls beyond 
some critical energy density s c , the droplet D hadronizes 
into an n-hadron configuration K = {h\h'2 ■ ■ - h n } with 
a probability proportional to £1, where SI represents the 
microcanonical partition function of an n-hadron system. 
Due to the huge configuration space, sophisticated meth- 
ods of statistical physics [22,23] have to be employed to 
solve the problem without further approximations. 

So our approach amounts to treating high density re- 
gions (droplets) for some time (between formation tj and 
hadronization r^) macroscopically, whereas before tj and 
after a microscopic treatment is employed. What hap- 
pens - microscopically - between tj and is not spe- 
cified, there may be a first or second order transition, just 
a crossover, or even some nonequilibrium transition. The 
macroscopic treatment is chosen due to the lack of ap- 



propriate transport theories of dense hadronic or quark 
matter. So at present we parametrize the behaviour of 
the dense matter in a simple fashion, the time evolution 
as longitudinal expansion and a hadronization according 
to n-body phase space. The hadronization of a droplet 
is not meant to represent a dynamical description of a 
phase transition, it means that at one observes a mul- 
tihadron system, whatever happened between tj and r^. 
Whether our parametrization is realistic, and what hap- 
pens microscopically between tj and ramains to be 
investigated by before mentioned theories. 

The first stage of our approach is the identification of 
high energy density regions, based on the string model, 
which is already discussed elsewhere [20]. Due to the 
empirically found correlation, y = £, between the aver- 
age rapidity y of particles and of the space-time rapid- 
ity C, a hypersurface % T of constant proper time r can 
be introduced, in the central region simply defined by 
t 2 -z 2 = t 2 . After having used the string model (VENUS 
5.08) to get complete information on hadron trajectories 
in space and time, we may now, for given r, determine 
energy densities on % T and thus locate high density re- 
gions on % T . 

High density regions are considered as QM droplets, 
presently it is assumed that they expand purely longit- 
udinally. Whenever other droplets or hadrons cross its 
way, the two objects fuse to form a new, more energetic 
droplet. Due to the expansion, the energy density of a 
droplet will at some stage drop below s c , which causes 
hadronization, to be described in the following sections. 

We consider the concept of QM droplets to be crucial, 
in particular at SPS energies. It has been shown [21] 
that at these energies energy density fluctuations are im- 
portant: one observes intermediate size regions of high 
density rather than a uniform distribution. Typical sizes 
of few tens of fm 3 are observed for these high density 
regions. 
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II. HADRONIZATION ACCORDING TO N-BODY 
PHASE SPACE 

For the hadronization of QM droplets we employ the 
following procedure: the probability of a droplet D with 
invariant mass E and volume V to hadronize into a con- 
figuration K = {hi, . . . , h n } of hadrons hi is given as 



prob(7> A) ~ fi(A') 



(1) 



with Sl(A') being the microcanonical partition func- 
tion of an ideal, relativistic gas of the n hadrons hi. 
We first have to define a set S of hadron species; we 
take S to contain the pseudoscalar and vector mesons 
(71", K, T], rj' , p, K* , ui, <f>) and the lowest spin--| and spin-| 
baryons (N, A, £, S, A, £*, S* , £1) and the corresponding 
antibaryons. A configuration is then an arbitrary set 
{hi, . . . ,h n ] with hi E S. 

The partition function is given as 



fi(A') — C V ol CdegCident 4> 



(2) 



with 



a 



V" 



vol 



(2TTh) 3 



n 1 
Cdeg = Y\. 9i ; Cident = J^J 1 



(3) 



Here, Cd eg accounts for degeneracies (gi is the degeneracy 
of particle i), and Cedent accounts for the occurrence of 
identical particles in K (n a is the number of particles of 
species a). The last factor 

<f> = <f)(E,mi, . . .,m n ) (4) 



is the so-called phase space integral, with Si = \/ mf + pf 
being the energy and pi the 3-momentum of particle i. 
The term ^q,s ? , ensures flavour conservation; qi is the 
flavour vector of hadron species i, and Q is the flavour 
vector of the droplet (the components of the flavour vec- 
tors represent the net quark content for the quark fla- 
vours u, d, . . .). The expression eq. (4) is valid for the 
centre-of-mass frame of the droplet D. 

We are going to employ Monte Carlo techniques, so we 
have to generate randomly configurations K according to 
the probability distribution Sl(A'). We want ot develop a 
method in particular for intermediate size droplets, cov- 
ering droplet masses from few GeV up to 100 or 1000 
GeV. So the method should work for particle numbers 
n = \K\ between 2 and 10 3 , which means, we have to 
deal with a huge configuration space. Such problems 
are well known in statistical physics, and the method at 
hand is to construct a Markov process, specified by an ini- 
tial configuration A'o, and a transition probability matrix 
p(Ki —> A'i_|_i). In generating a sequence A'o, Ki, K2, ■ ■ 
two fundamental issues have to be payed attention at: 



• initial transient: starting usually off equilibrium, it 
takes a num 
equilibrium 



takes a number of iterations, 7 eq , before one reaches 



• autocorrelation in equilibrium: even in equilibrium, 
subsequent configurations, K a and K a +i, are cor- 
related for some range /auto of i. 

In general, both 7 eq and J a uto should be as small as pos- 
sible. 

We are going to proceed as follows: for a given 
droplet D with mass E and volume V , we start from 
some initial configuration A'o, and generate a sequence 
A'o, A'i , . . . , A'/ e , with 7 eq being sufficiently large to have 
reached equilibrium (which is defined to be the steady 
state of the Markov process). If we repeat this procedure 

many times, getting configurations A'j 1 - 1 , A'j 2 - 1 , . . . , these 
configurations are distributed as Sl(A'). So for our prob- 
lem, we have only to deal with the initial transient, not 
with the autocorrelation in equilibrium. We have to find 
a transition probability p such that it leads to an equi- 
librium distribution Sl(A'), with the initial transient 7 eq 
being as small as possible. 

So our task is twofold: we need to find efficient ways 
to calculate, for given A', the partition function Sl(A'), 
and we have to find an appropriate transition probability 
p(K a —> K b ). 



III. THE PARTITION FUNCTION Q(K) 
The partition function is given as (eq. (2)) 

fi(A') = CVol Cdeg C ld ent <j> (5) 

with the phase-space integral <f> (eq. (4)) and some pre- 
factors C'i (eq. (3)). In the following we discuss methods 
to calculate <f> for an arbitrary number n of particles, 
starting with n = 2. 

For 2 particles (n = 2), we have 



<f>{E, mi , m 2 ) = j d 3 pi d 3 pi 8{E - Si - s 2 ) 6(p 1 + p 2 ) 
= 47r / dpp 2 8(E — \Jm\+ p 2 — \jm\ +p 2 ). (6) 



With po representing the root of the argument of the 
(^-function, 



Po = 
we get 



£2_ 2(m 2 +m 2 ) + i_ (m 2_ m 2 )2 



1/2 



(7) 



<f){E,m 1 ,m 2 ) = A'Kp Umf + p A ) = + (mi, + p A ) = 



(8) 
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In the following we consider n > 3. We propose a 
method to calculate <f>, introduced by Hagedorn [24]. The 
phase-space integral, eq. (4), may be written as 

<j>(E, mi, ... , rn n ) (9) 

{ . n n n 

= (4^) n j n n s ( e - e ^) w ^ .•••.?"). 

i = 1 z = 1 i = 1 

with pi = | pi | , and with the "random walk function" W 
given as 

i /■ n n 
W ^---,Vn):=j^ JJldeiSi^piei), (10) 



with i{ = Pi/\pi\, and with 

dii = sin $i ddi dipi 



(11) 



representing the integration over all directions for a given 
length \ fi | of a momentum vector of a particle. The name 
"random walk function" is due to the fact that W rep- 
resents the probability to return back to the origin after 
n "random walks" pi 'ti with given step sizes pi. 
We first evaluate W for n = 3. One may write 



3 ( 3 

w( P1 ,p 2 ,P3)= f n^?* jn 



-4^2- f*( S «0. ( 12 ) 



with cii := The integration over d 3 Q3 may be per- 

formed, and one obtains 



W(pi,p 2 ,P3) 



1 1 



(13) 



(4tt) 3 pXpIpI 

x y rf 3 gi rf 3 g 2 K a \ -Pi) H a 2 -P'z) <5(l?i + q-z\ ~ Pz)- 
Taking $ to be the angle between qi and q2, we have 

|Si + 92I = \/ a i+ a 'i~ 2aia 2 cos 

= \/p'{+pI- 'IpiP'i cos «?, (14) 
so we may perform five of the six integrations, to obtain 
2(2tt) 2 1 



W(pi,p 2 ,P3) 



(4tt)3 pi 



(15) 



and we get the final result 

W( Pl ,P2, P 3) (16) 
(87T P1P2P3)" 1 if |Pl - P2I < P3 < \Pl + P2| 







else 



In the following, we discuss methods to calculate W, 
for n > 4. The random walk function may be written as 



W(pi,...,Pn) 
which leads to 

W(pi ,...,p r 



1 1 

(4tt)" (27r) a 



a / 



7 C 1 



■1 /* CO n 



sin A 

1 Pi X 



(17) 



(18) 



This is easy to evaluate via numerical integration, as long 
as n e ff is large (> 10), where n e fj is the number of mo- 
menta with pi > e, with some small e. Otherwise the 
integrand fluctuates so strongly, we have to find a differ- 
ent method, as discussed in the following. 
From eq. (18), we get 



W( Pl ,..., Pn ) = — 



4tt 2 U Pj (2i) n 
00 ~ i£ dX 



x 



A n-2 



(19) 

n 



The product Y[{} m ec l- (19) can be written as 



n e ^ ixasps , 

j = l a^il-l} 



(20) 



which is equal to 



J2 <ri ■ ■ ■ <r„ e iAS<7 ^' . (21) 

o\ . ..a n 



For ^(TjPj being nonnegative, we have 

noo — ie 



d\ 



X n-2 

= 2tt» Res{A 2 - n e iAS<7 ^} A=0 



1 f rl"~ 3 

2«— 4^e* xs ^Pi 



(n-3)! [dX 



A = 



(22) 



whereas for negative Y,(jjPj, the integral is zero. So we 
get 



x / d cos?? 6(\Jpl + P2 ~ 2piP2 costf -p 3 ), W(pi , p n ) 



1 



2 n +! 7r(n - 3)!pi...p„ 

X ^ (7i...<7„(^( 



\ n — 3 



(23) 



To perform the summation J2 CT CT , it is useful to 
take a specific sequence of it's (using a = {<J\, ■ ■ ■ , & n }), 
namely [24] 
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?(i) 

?(2) 
?(3) 
f(4) 



{+ + + ■ 

{- + + 

{-- + 
{ + 



•} 
•} 

+ +...} 
+ +...} 



(24) 



where we simply write + and — rather than +1 and —1. 
The general rule is that (f( v+1 ) is obtained from tj( v \ by 
changing position number j v , where 31,32,33 • • • is given 
as 



1,2,1,3,1,2,1,4,1,2,1,. 



(25) 



with obvious continuation; the rule to obtain the se- 
quence {j v } is: taking the binary representation off, one 
obtains j v as the position of the right-most non-zero digit, 
counting from right to left. For example for v = 7 = 111, 
we have jV = 1, for v = 8 = 1000, we have j$ = 4. With 
this prescription actually all possible it's are accounted 
for. Using this sequence of it's as described above, we 
obtain 



J2 a i p i 



J2 a i p i 



2^lW, (26) 



<j(") 



which means that rather than performing the whole sum 
^2<TiPi, only one term, —2a^pj u , has to be added. 

Eq. (23) provides a method to calculate W, as long 
as the number n of particles is not too large (n < 20). 
We discussed earlier (eq. (18)) a way to calculate W for 
large n (n e ff > 10). Fortunately, there is some overlap 
between the two methods, so we may always use one or 
the other procedure to evaluate W to any given accuracy. 
In practice, we use eq. (23) for n < no (with no — 10), 
for larger n we use eq. (18). It may happen that for 
n > no the desired accuracy cannot be achieved, due to 
the fact that one or several momenta p,' are small leading 
to a strong oscillation of the integrand of eq. (18). In 
this case, we use the other method, eq. (23). 

Having a reliable and efficient method to calculate W, 
we may return to the problem of how to calculate the 
phase space integral <f> efficiently. From eq. (9), we obtain 

/'OO /'CO n 

4>(E, mi, ... , m n ) = (47r) n / dei . . . de n J^J pi a 

n 

x6(E-J2ti)W(p 1 ,...,p n ), (27) 

8 = 1 

af ter a cha nge of variables toward particle energies e% = 
\J mf + pi . We introduce kinetic energy variables, 



ti . — Si fn^ , 
and a total kinetic energy T, 



(28) 



(29) 



and obtain 

poo poo n 

<f>(E,m 1 ,...,m n ) = (ATr) n dtx . . . \ dt n TT p { 

Jo Jo i=1 

n 

x6(T-J2ti)W(p 1 ,...,p n ). (30) 

8 = 1 

We now introduce "accumulated" kinetic energies s 8 - via 



E 



(31) 



with the inverse 

ti = Si - s 8 '_i, s = 0. (32) 
The variables are replaced successively, 



t\ = s\ , dt\ = ds\ , 

t 2 = s 2 - si, dt 2 = ds 2 , 

tn — &n Sn — 1, dt n — d-S n , 



(33) 



and we obtain 

poo poo poo 

<f)(E,mi, . . .,m n ) = (47r) n / dsi ds 2 . . . / ds n 

JO J S 1 J S n - 1 



Y[p i e i 6(T-s n )W(p 1 ,...,p n ). (34) 



The integration over s n is trivial and may be performed, 
to obtain 

<f)(E, mi, ... , m n ) 

poo poo poo pT 

= (47r) n / dsi ds 2 . . . ds n _ 2 / cten-i 

JO J Si J S n -3 J S n -2 



Y[pi£i W(pi,. ..,p n ). 



(35) 



All upper limits may be replaced by T. Introducing en- 
ergy fractions, 



— £i 

X i '■— rp J 



(36) 



we get 



^(E,m 1 ,...,m n ) = (^) n T n - 1 (37) 

n 

dxi . . .dx n -iY\.Pi £ i W (Pl> ■ ■ - ,Pn)- 



'0<ii<..Xi,_i<l 

Using the definition 

/^jj-Xn rpn-1 ™_ 
1p(Pl>--->Pn) ■= , _ l y Y[pi£iW(pi,...,p n ), 

\ n >' 8 = 1 



(38) 
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we may write 
<j>(E, mi, ... , m n ) 

= (n-iy. 



(39) 

dxi . . .dx n -iip(xi, . . ., aj„_i), 



0<ii<...<i„_i<l 

where ip(xi , . . . , £ n -i) is meant to be tp(pi , . . . , p n ) with 
Pi and a expressed in terms of x\, . . .,x n _\. This may 
be solved via Monte Carlo as 

1 N 

4{E, mi,..,m„) = -^ rP(x[ P) . . . x^l,), (40) 



13 = 1 



where the x^ are ordered random numbers, 



< < 4^ < . - . < < 1- 



(41) 



So for each Monte Carlo step, n — 1 random numbers have 
to be generated, ordered according to size, and then used 
to evaluate ip(x l 'f\...,x l ' 1 j 3 } 1 ). To avoid ordering, one 
may introduce the variables 



(42) 



using the definition x n := 1. We get 

n-1 

j=i+l 

the last equation holding for i < n — 1; so we have 

n — 1 n — 1 n — 2 n — 1 



(43) 



n = n dzi n n 

z' = l z' = l i = l j=i-\-l 

n — 1 n — 1 



(44) 



From eq. (39), we get 
<f)(E,mi, . . .,m n ) (45) 

i-l i-l "-1 

= / dzi . . . J dz n _ 1 T\iz i *~ 1 ip(z 1 ,...,z n _ 1 ), 
■Jo Jo „•_ -, 



where obviously t/"( z i ; • • • ; z n-i) is meant to be 
ip{pi, ■ ■ ■ ,Pn) with pi expressed in terms of Z{. We now 
introduce 



n := / i?- i d£ = z i t 
lo 



(46) 



to obta 



4>(E, mi, ... , m n ) = I dn... dr n _i t/>(ri , . . . , r n _i). 



o 



o 



(47) 



The ri are now uncorrelated, no ordering is required. A 
Monte Carlo solution is simply 



SP) 



1 N 

<f > (E,m 1 ,...,m n ) = -J2Hr ( f\---,ri P > 1 ), (48) 

13 = 1 

with uncorrelated random number rf \ So for each 
Monte Carlo step [3 the following procedure is followed 
(we drop the index f3): 

• generate n — 1 random number r 8 -; 

• calculate z 8 - = tfrl and then the energy fractions 

Xi — Xi^.\Zi, X n — 1, 

• calculate the accumulated energies s 8 - = Txi, and 
then the kinetic energies ti = s 8 — s 8 _i, using s = 0. 
Then calculate t he energies g = ti + rrii and the 
momenta pi = \/ti(ti + 2m,-); 

• calculate tp(pi p n ) according to eq. (38), by us- 
ing the above methods to calculate W{p\, . . ■ ,p n )- 

Summing up all the t/>'s and dividing by the number 
N of Monte Carlo iterations (see eq. (48) provides 
the Monte Carlo result for <f>(E, mi, m n ). Clearly 
most of the computing time goes into the calculation 
of W(pi , . . . , p n ), in particular into the evaluation of 
sinpiX/piX for large n (see eq. (18)). 



IV. THE METROPOLIS ALGORITHM 

As mentioned earlier, we want to generate randomly 
hadron configurations K = {hi, . . . , h n } according to the 
probability distribution Sl(A'), where Sl(A') is the mi- 
crocanonical partition function discussed extensively in 
the previous section. With A'( a ) being such configura- 
tions, mean values of observables O(A') are then simply 
calculated as 



< O >-- 



1 N 



(a) 



(49) 



To construct a configuration K^ a \ for each a, a chain of 
configurations A'o, A'i, A'2, . . . , A'/ e is constructed, which 
is characterized by an initial configuration A'o and a ran- 
dom matrix p( A'i — ► A' 8 _|_i), which specifies the probab- 
ility of a configuration A' s - being followed by A' 8 _|_i. The 
number I eq of iterations must be large enough to ensure 
equilibrium, only then the randomly generated A'/ eq are 
distributed as Sl(A'), for an appropriate p, and we take 



A'(«) = A, 



(50) 



In the following, we discuss how to construct an "ap- 
propriate" p, which makes the A'( a ) being distributed as 
fi(A'). 
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Sufficient for the convergence to Sl(A') is the detailed 
balance condition, 



Q(K a )p(K a -+ K b ) = Q(K b )p(K b -+ K a ) 



(51) 



and ergodicity, which means that for any K a , K b there 
must exist some r with the probability to get from K a 
to Kb in r steps being nonzero. Henceforth, we use the 
abbreviations 

Cl a := tt(K a ); pab := p(K a —> K b ). (52) 

Following Metropolis [22], we make the ansatz 

Pab = W ah U a b , (53) 

with a so-called proposal matrix w and an acceptance 
matrix u. Detailed balance now reads 



Ugb _ &b W h g 
U ha &a W a b 

which is obviously fulfilled for 



(54) 



Uab = F 



&a Wab 



(55) 



with some function F fulfilling F(z) / F(z 1 ) = z. 
lowing Metropolis [22], we take 

F(z) = min(z, 1) . 



Fol- 



(56) 



The power of the method is due to the fact that an arbit- 
rary w may be chosen, in connection with u being given 
by eq. (55). So the task is twofold: one needs an efficient 
algorithm to calculate, for given K, the weight Sl(A'), 
and one needs to find an appropriate proposal matrix w 
which leads to fast convergence (small I eq ). The first task 
can be solved, as shown in the previous section. In the 
following we discuss about constructing an appropriate 
matrix w. 

Most natural, though not necessary, is to consider sym- 
metric proposal matrices, w a b = Wb a , which simplifies the 
acceptance matrix to u a b = F(£li,/Cl a ). This is usually 
referred to as Metropolis algorithm. Whereas for spin 
system, it is obvious how to define a symmetric matrix 
w, this is not so clear in our case. We may take spin 
systems as guidance. A configuration K is per def. a set 
of hadrons {hi, . . . , h n } with the ordering not being rel- 
evant, so {ir° , 7T° ,p} is the same as }. We intro- 
duce "microconfigurations" to be sequences {hi, . . . , h n } 
of hadrons, where the ordering does matter. So for a 
given configuration K a = {hi, . . . , h n } there exist several 



microconfigurations K 



{hjr. 



(i)> 



}, with 7Tj 



representing a permutation. The weight of a microcon- 
figuration is 



Cl(Ka) 



(57) 



with n a being the number of hadrons of type a. Taking 
for example K = {p,ir° ,ir }, there are three microcon- 
figurations {p, 7r°,7r }, {ir°,p, it } and {7r°,7r°,p}, with 
weight fi(A')/3. 

So far we deal with sequences {hi, . . . , h n } of arbitrary 
length n, to be compared with spin system with fixed lat- 
tice size. We therefore introduce "zeros", i.e. we supple- 
ment the sequences {hi, . . . , h n } by adding L — n zeros, as 



{hi 



h n ,0, 



, 0}, to obtain sequences of fixed length 



L. The zeros may be inserted at any place, not necessar- 
ily at the end. Therefore the weight of a microconfigur- 
ation K a j with zeros relative to the one without, K a j, is 
one divided by the number of possibilities to insert L — n 
zeros, so from eq. (57) we get 



1 



n 

.a£S 



n\ (L — n)\ 



Cl(Ka) . (58) 



We now have the analogy with a spin system: we have a 
one-dimensional lattice of fixed size L, with each lattice 
site containing either a hadron or a zero. Henceforth, 
we use for microconfigurations with zeros the notation 
K a j = {hi, . . . , /il} with hi being a hadron or zero. 

Since from now on we only consider microconfigura- 
tions with zeros (K a j) rather than configurations (K a ), 
we are going to write K a instead of K a j , keeping in mind 
that a represents a double index, and say "configuration" 
rather than "microconfiguration with zeros" . The ad- 
vantage is that we can use the above formulas specifying 
the Metropolis algorithm without changes. 

We are now in a position to define a symmetric pro- 
posal matrix w(K a —> Kb), with K a = {h", . . . , h a L { and 
K b = {h\,...,h b L }, as 



w(K a —> K b ) 



(59) 



L(L-l) 



with 



v(hw ^ Kh\) 



\Vihfh])]- 1 if h\h) e V(hfh]) 
else 



(60) 



where V(hfhj) is the set of all pairs (hihj) with the same 
total flavour as the pair (h"hj). The symbol \V\ refers 
to the number of pairs of V. The term {} in eq. (59) 
makes sure that up to one pair all hadrons in K a and 
Kb are the same, the term L(L — l)/2 is the probability 
to randomly choose some pair of lattice indices i and j. 
So our proposal matrix amounts to randomly choosing a 
pair in K a , and replacing this pair by some pair with the 
same flavour, with all possible replacements having the 
same weight. The proposal matrix is obviously symmet- 
ric, since v is symmetric (the symmetry of v is crucial!). 
We have now fully defined an algorithm, which due to 
general theorems will converge, but how fast, i.e., how 
large is 7 eq ? This is going to be investigated later. 
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V. A GENERALIZED CONFIGURATION SPACE 

So far, a configuration was defined to be given as 
K = {hi, . . . ,h n } with hi specifying hadrons species, 
for example K = {7r°,7r°,p}. The Metropolis Algorithm 
introduced in the previous section will provide random 
configuration K^ a \ which are distributed as Sl(A'). This 
approach is not yet satisfactory for the following reasons: 
we do not want to make predictions for multiplicities 
only, but also consider momentum distributions of the 
hadrons; the method is also extremely slow due to the 
fact that, for each Metropolis step, the function Sl(A') 
has to be evaluated, which itself requires a Monte Carlo 
procedure with many iterations. There is a way to cure 
both problems: one has to consider a generalized con- 
figuration space, such that not only hadron species are 
considered but also hadron momenta. 

A naive generalization would be to introduce config- 
urations as {hi, . . . , h n ;pi, . . . ,p n }, with pi representing 
the particle momenta. The symbols hi represent again 
the hadron species. There are two problems about the 
naive generalization: the momenta are not independ- 
ent, since their sum must be zero, and, in addition, 
we loose the symmetry property of the proposal matrix 
w(K a —> Kb)- This symmetry is not really necessary, but 
at least one needs to be able to calculate the asymmetry 
w(K a —> K h )/w(K h —> K a ). 

To find a reasonable generalization one should recall 
the discussion following eq. (27), where a couple of co- 
ordinate transformations were applied to calculate the 
phase space integral <f>. The final result was eq. (47), 

<j>(E, mi, ... , m n ) (61) 

dri . . . / rfr„_i ip(E, m 1 , . . . , m n ; 7*1, . . . , r n _i), 
/o Jo 

with ip given in eq. (38) as 

ip(E, mi, . . . , m n ; ri, . . . , r n _i) (62) 

( 47r )n T n-l 



where the r- are related to the momenta p,' via eq. (63). 
The weight of such a configuration is 



YrpII K£i w(pi,...,p n ). 



Here we also indicate the dependence of ip on E and 
mi,...,m n , which has been dropped in the previous 
chapter. The symbol T denotes the total kinetic energy 
E — Em;, and the absolute values of the momenta are 
expressed in terms of the r- as 



Pi = \/t{(ti + 2m-) 

ti = T(xi - Xi_i), x = 

Xi — Xi^-l^Jv~i, X n — 1. 



(63) 



Contrary to the pi, the r- are independent of each other. 
Based on eq. (61), we introduce generalized configura- 
tions G as 



11(G) — C V ol Cdeg Cident V" 



(65) 



(see eq. (2)), with ip given in eq. (62). We always use the 
same symbol £1 for the different functions Cl(x), depend- 
ing whether a; is a configuration K a , a microconfiguration 
K a j or a generalized configuration Ga- 
in order to define a symmetric proposal matrix 
w(K a —> Ki,), we introduced in the previous chapter 
microconfigurations. We proceed similarly for the 
generalized configurations G a . For a given G a = 
{hi, h n ; ri, r n _i}, one obtains several microconfig- 
urations G a j, by introducing L — n zeros, leading to a 
sequence of hadrons {hi, h^} of fixed length, with hi 
being a hadron or zero. The sequence ri, r n _i repres- 
ent the momenta of the nonzero hi (the hadrons). We 
supplement ri, r n _i by L — n numbers r n , . . . , rt_i, 
with < r- < 1. A generalized microconfiguration is 
thus given as 

G a j = {hi, . . . , h L ; ri, . . . , r L _i}, (66) 

with weight 

^(C a j) = C vo l Cdeg Cident C m ; cro if), (67) 

with ip = il>(E , mi, . . . , m n ; ri, . . . , r n -i) given in eq. 
(62), the prefactors C vo i, C deg , C lden t given in eq. (3), 
and the other prefactor given as 



Cry 



n\(L — n)\ 




(68) 



G = {hi 



<h n ;ri, 



i}, 



(64) 



(see eq. (55)). The r- for i > n seem to be obsolete, since 
only 7*i, . . . , r n _i are needed to calculate the momenta of 
the n hadrons, however, the numbers are needed to define 
a symmetric proposal matrix. As for the configuration K , 
also for the generalized ones, we write simply G a rather 
than G a j and drop the term "micro" . 

We are now going to define a symmetric proposal mat- 
rix w(G a —> Gb), with G a given as 

G a = {h a i,...,h a L ;r a 1 ,...,rl_ 1 }, (69) 

and Gj correspondingly. We introduce a "species part" , 

K a = {K,...,h a L }, (70) 

and a "momentum part" , 

R a = {r a i,...,r a L _ 1 }, (71) 

and corresponding definitions for A'j and We may 
now define a proposal matrix w as 

w(G a — > Gt) = w spec (K a — > Kb) 

•fmom (Ra^Rb), (72) 

where the "species matrix" u> spe c is defined in eq. (16), 
with w instead of u> spe c being used. The "momentum 
matrix" is defined as 
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. iVmom-l 



(R a —> Rb) — Y\. ^Rc, W mom(Ra ~> R Cl ) 

■' 8 = 1 

X "'momf^ci ~^ Rc 2 ) • • • 

x ■■•in,(V„-ii ^Rb), (73) 



with dR = dridr 2 ■ ■ .dri-i, and with 



1 L ~ 1 



4). 



(74) 



The term (L — 1) _1 indicates the probability to randomly 
choose a position i between 1 and L — 1, the second part 
of eq. (74) ensures that all r? for j ^ i are not allowed 
to be changed, however, the number rf is replaced by an 
arbitrary number rf £ [0, 1] with probability one. The 
following Monte Carlo procedure generates an "updated" 
Rb, starting from R a , according to eq. (74): choose ran- 
domly a position i(l < % ' < L — 1), and then replace rf by 
some random number r £ [0,1]. This provides Eq. 
(73) simply accounts for repeating the above procedure 
ffmom times. 

Nmom is an important technical parameter of our pro- 
cedure (in addition to J eq and L), which may be chosen 
between 1 and L — l. Let us consider the weight £1(G) for 
fixed hadron species hi, . . . , hi, but varying momentum 
variables R = {ri, . . . , rt_i}. Out of the huge R phase 
space (for large n) only a very small region contributes 
with significant weight; for most values of R, £1(G) is 
practically zero. So taking N moTn = L — 1, representing a 
complete i?-update, would frequently propose configura- 
tions with zero weight, which are rejected with a large 
probability. So, one may get trapped for a long time. 
Clearly, this choice of N moTn leads to large equilibration 
times J eq . The other extreme, N moTn = 1, provides up- 
dated configuration very close to the original one. Now 
it takes a long time to test the available phase space. In 
particular it might easily happen, that one gets trapped 
in the neighbourhood as a local maximum. We have ac- 
tually the following situation: for a given number n of 
hadrons in G, we have a local maximum of £1(G) at some 
G^ n \ with £1 dropping very fast with G moving away from 
G in \ The maximum values fi(G W ) for 

neighbouring 

n's are not so different though. One easily gets trapped 
around some G 1 -"- 1 , even with n being quite far away from 
the equilibrium value. So N moTn must be chosen large 
enough to explore the available phase space without get- 
ting trapped at a local maximum, but not too large, to 
avoid exploring extremely unlikely regions. 



VI. AN ASYMMETRIC PROPOSAL MATRIX 



Considering particle ratios, like n w o /n w +, we find im- 
mediately that we have a very slow convergence, so 7 eq is 
too large for the method to be of practical importance. 
This is obvious, since the proposal matrix w does not 



Kv 



h b 


h b 


h b 

— f 


m 


n 





FIG. 1. Double pair exchange. 

act very democraticly: flavourless particles like 7r° , p° or 
also zeros are much more frequently proposed than all 
the rest. This shortcoming can be fixed by defining w 
such that two pairs are exchanged rather than one, the 
first pair being replaced by a completely arbitrary pair, 
the second one by some pair to guarantee flavour conser- 
vation. 

Such a proposal matrix is not symmetric any more, and 
the new method is therefore referred to as "asymmetric" 
or "double-pair-exchange" procedure. In the following, 
we provide some details about the asymmetric method. 

The new proposal matrix is still of the form w = 
K'spec'fmom, see eq. (72), where u> spe c refers to particle 
species and w mom to momenta. We take the same w mom 
as before, only u> spe c is changed. We define 



w spec (K a —> K h ) 



4 ; (i + i5i 



£ n ^ 

m<n<z'< < ; I k = 1 



xvdhlhlhirhlhfh^^h'h'), (75) 
with h representing the antiparticle of h, and with 

v({hi,h 2 ,..]^h\h]) (76) 

C if V({hi,h 2 ,...}) empty 

= ^ \V({hi,h 2 ,...})\- 1 iih\h)£V({hi,h 2 ,...}) . 
{ else 

V({hi, h 2 , ■ ■ ■}) represents the set of all pairs h{hj of 
hadrons with the same flavour as the set of hadrons 
{hi, h 2 , . . .}. So the double pair exchange works as fol- 
lows (see fig. 1): two pairs m < n and i < j (with n < i) 
are chosen randomly, with equal probability ( 4 ) for all 
possible double pairs. The first pair h a n h a m is replaced by 
some arbitrary pair h h n h h m , all possible pairs having equal 
weight (1 + |iS|) -2 , with |<S| being the number of had- 
rons in the basic hadron set S (containing the standard 
hadrons, but for testing purposes we will later also use 
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£3/4 : 16 



J_ 

16 



Kv 



FIG. 2. Example for double pair exchange. 

reduced hadron sets). We have 1 + |<S| rather than |<S|, 
since we are also considering "zero" . To ensure flavour 
conservation, the second pair hfh" has to be replaced by 
some pair h\h h j having the same flavour as the set of had- 
rons {h^, h®, h h m , h h n , hf , hj}, where all possible pairs are 
taken with equal weight. 

The new proposal matrix is no more symmetric, which 
means, the acceptance matrix has to be calculated ac- 
cording to the general expression 



u(G a —> Gt) = F 



w(G b ^ G a ) Cl(G b ) 
w(G a -+ G b ) fi(G a ) 



(77) 



rather than using simply u(G a —> G b ) = 
F({l(Gi,)/{l(G a ))- Since u> mom is symmetric, the "asym- 
metry" is given as 



w(G b — » G a ) _ w spec (Kb — > Kg) 
w(G a —> G h ) w spec (K a —> K h ) ' 



(78) 



To evaluate the r.h.s. of eq. (78) we simply need to 
calculate the ratio of the probabilities v to exchange the 
second pair, which is given as 



\v({hi,hih° m ,h a n ,hihn)\-i' 



(79) 



Let us discuss a simple example (see fig. 2). After 
choosing randomly positions m < n < i < j, we may 
find for example two pairs each consisting of two 7r°. So 



we nave 



h a 



h a 



7r . The first pair 



is replaced by some arbitrary hadron pair, each pair is 
considered with equal probability (l + |iS|) -2 . Lets choose 
a (tt + , it ) pair. In order to achieve flavour conservation, 
the new second pair must have the same flavour as the 
set of hadrons {h^, hf,h h m ,h h n , hf , hj} , which is in this 
case {it , it , tt + , 7r°, 7r°, 7r }, so the flavour is ud. Taking 
just for testing purposes a reduced hadron set 1S3/4 = 
{it , 7r + , 7t - }, the set of pairs with flavour ud is 



{(0,7r-),(7r°,7r-),(7r-,0),(7r-,7r )}. 



(80) 



Taking equal probabilities, the weight to choose any of 
these, for example (7r _ ,7r°), is 1/4. Taking the inverse 
case, we have the two pairs (h h m ,h h n ) = (7r + ,7r°) and 
(hf,h b j) = (7T - ,7r°) where the first pair, (7r + ,7r°) is re- 
placed by (hf n) hf l ) = (7r°,7r°), with the probability for 
this replacement being (1 + |i5|) -2 = 1/16. The second 
pair, (7T - ,7r°), has to be replaced by (hf,hj) = (7r°,7r°), 
but the probability for this is not 1/4. How many 
pairs would be possible? The pair must have the fla- 
vour of the set {h h m , h^jh^jh®, h\ , h b A , which is here 



{tT 1 



7T , 7T , 7T , 7T 



V}, 



the flavour must be 0. The 



set of possible pairs is 
{(0, 0), (0, 7T ), (tt , 0), (tt , tt ), (tt+tt-), (81) 

so the probability to select a pair is 1/6. The asymmetry 
w(G h -+ G a )/w(G a -+ G b ) is therefore (l/6)/(l/4) = 
2/3. 

The example demonstrates that, indeed, the proposal 
matrix is in general asymmetric, however the asymmetry 
can be calculated quite easily. For counting the number 
of possible pairs, one just has to make sure to account for 
the ordering, for example (7r + , 7r~) and (tt - , 7r + ) must be 
considered as different pairs. 

The basic set S of hadrons has been defined to con- 
tain mesons and (anti)baryons from the two lowest mul- 
tiplets each. For testing purposes we introduce "test set- 
s'' iS^, for example, we define S\ = {"Tr }, with the 7r° 
being considered massless, we use S2 = {i" } as well, 
but considering the correct mass. The complete list 
iSi,c>2, • • - ,S9,Sio = S is given in table I. We consider 
test sets with massless hadrons, because in this case an 

TABLE I. The hadron sets <S M . Odd fi implies massless 
hadrons, for even fi the correct masses are considered. 



1,2 
3,4 
5,6 
7,8 
9,10 



hadrons in <S M 
IP 

— + — 

7T = 7T ,7T ,7T 

ir,p,n,p,n 

7r, if ,»7,»/,|),ra,E,H & antibaryons 

as 7,8 & p, K* , lo,4>,A, £*, £*, «~ & antib. 



analytical treatment is possible, providing useful checks 
of our Monte Carlo procedures. We will discuss the ana- 
lytical treatment and detailed comparisons between ana- 
lytical and Monte Carlo results later. Presently, we are 
just interested how fast our asymmetric algorithm leads 
to convergence, depending on the size of the hadron set. 
We consider a "test droplet" of size V = 10 fm 3 with 
mass E = 10 GeV, and we apply our hadronization pro- 
cedure for the test set Si , 1S3, 1S5, 1S7, and 1S9 (so we restrict 
ourselves to massless hadrons). In fig. 3 we plot the mul- 
tiplicity n versus the number of iterations, for different 
sets Sfi . One clearly observes a fast convergence for small 
sets, but for 1S7 and in particular 1S9 (containing the full 
set of hadrons, just massless), we have a very slow con- 
vergence. We discuss in the next section a method to 
improve that. 
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25 



20 



15 



10 



E = 10.0 V = 10.0 




1000 2000 3000 4000 
iterations 



FIG. 3. The multiplicity n versus the number of itera- 
tions, for different sets 5 M : 5i (solid line), 53 (dotted), 57 
(dashed), 59 (dash-dotted). An asymmetric proposal matrix 
is employed. 



VII. A VERY ASYMMETRIC PROPOSAL 
MATRIX 

So far we have a method which converges fast for 
small hadron sets (Si , S3) but very slowly for the large 
ones [St,Sq), and unfortunately the largest is the real- 
istic case. How can one improve the method? We re- 



call that "0" is treated like a hadron: when proposing 
a new pair hi,h,2, each hi may be any hadron from S^ 
or "0", so one may consider extended sets S^, contain- 
ing the hadrons from S^ and in addition "0" . So we 
have S\ = {0,7r }, S3 = {0, 7r°, 7r + , 7r - } and so on. A 
major difference between different S^'s is, that the rel- 
ative weight of the zeros in S^ decreases with increasing 
\i: where as for Si the zero has weight 1/2, for Sg the 
weight is only 1/55. On the other hand, large weight 
implies a large probability to propose a pair containing 
one or even two zeros, which may reduce the multiplicity, 
whereas small weight for zeros implies a large probabil- 
ity to propose pairs without zeros, making a reduction 
of multiplicity rare. So for Sg or S7, with small weights 
for the zeros, there is a large asymmetry in exploring the 
phase space, it is much more likely to propose a configur- 
ation with increased multiplicity than one with reduced 
multiplicity. Such an asymmetry, causing many unsuc- 
cessful suggestions, leads to a slow convergence. 

It is obvious how to improve our method: in case of 
large sets S^, the weight of the "0" must be increased 
relative to the hadrons. Since in this way we introduce 
another asymmetry to the proposal matrix w, we refer 
to w as the "very asymmetric" proposal matrix, to dis- 
tinguish from the "symmetric" case (eq. (59)) and the 
"asymmetric" case (eq. (75)). The "very asymmetric" 
proposal matrix w is defined in the following. We again 
use w = iDspecWmom, see eq. (72), where u> mom is taken 
as before (eq. (73)), only the "species matrix" u> sp ec is 
changed. We define 



w spec (K a —> K b ) 



(82) 



1 



4; (# Z ero + \S\f 



e n *>v 



m<n<z'<i I k = 1 



xv{{h a m ,h a n ,h b m ,h b n ,K,h]}^h?h] 

with 

v{{hi,h 2) ..]^h\h]) = 



(83) 



if V({hi,h 2 ,...}) empty 



\\n{hi,h2,..})\\ 

else 



with 



Z(h) 



N, 



zero 

1 



if h 

else 



(84) 



V({hi, h'2, ■ ■ •}) represents the set of all pairs hihj of 
hadrons (or "0") with the same flavour as the set 
{hi, h'2, ■ ■ ■}■ The symbol \\V\\ represents a weighted sum 
of pairs of V , 



wv\\= J2 m)z{hi). 



(85) 
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5 - 
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asymmetry is given as 

mi^h^hi'hihf^mi- 1 



\\n\h 



b Lb La La fob 



(86) 



just with || ... || instead of | . . . |. The "zero weight" N zeTO 
(meant to be an integer larger or equal to one) is a tech- 
nical parameter, which has to be chosen to guarantee fast 
convergence for a given hadron set S, in particular for the 
fall set S9 (or <Sio). 

From the fact that the "asymmetric method" (which 
corresponds to N zeTO = 1) works well for Si, where the 
weight of the zero is 1/2, one might expect in general 
good results for N zeTO = liS^I, i.e. N zeTO = 54 for S9 
or e>io- In fig. 4 we see that, indeed, the performance 
improves significantly by increasing N zeTO from 1 to 54. 

We did finally set up an algorithm, which seems to be 
sufficiently fast also for a realistic hadron set (S9 and 
Sio). In the following sections we are going to compare 
the Monte Carlo results with analytical calculations. 



VIII. THE ZERO-MASS LIMIT 

The hadron masses crucially affect the actual results 
of the simulations. However, just in order to test the nu- 
merical procedures, it is useful to consider the "zero-mass 
limit" , i.e. the case of all hadron masses set equal to zero. 
In this case analytical results can be obtained, which may 
be compared with our Monte Carlo simulations. We in- 
troduced already, for testing purposes, several basic had- 
ron sets Si (see table I), where Si, S3, S5, S7, and S9 refer 
to massless hadrons, with an increasing number of had- 
rons considered. In the following we discuss the analyt- 
ical treatment for the case of massless hadrons [25,26]. 

We consider a configuration with the flavour of the n— 
hadron system, X2I=i lit being equal to the flavour Q of 
the droplet D. The phase space integral is then 



8 = 1 



FIG. 4. The multiplicity n versus the number of iterations, for Representing the 6-function as 
different values of N zevo : 1 (solid line), 3 (dotted), 54 (dashed). ^ 
We refer to a droplet with E = 10 GeV, V = 10 fm 3 and the 6(x) = 

(Z7T) 



hadron set 59. 

Whereas in eq. (75) all pairs h\h h j have the same relative 
weight 1, we now consider a relative weight Z(hf)Z(h b j), 
which means, the zero has a weight being N zeTO times 
larger than that of a hadron. Having defined the pro- 
posal matrix w, we have to determine the asymmetry 



and the ^-function as 



(x) = 

V ' 2th 



d 3 X e ixs 



da 



(87) 



(88) 



(89) 



the phase space integral may be written as 



1 



W spec (Kb —> Ka)/w spec (Ka - 

late the acceptance matrix. 



> Kb), necessary to calcu- 
Similar to eq. (79), the 



z"(2tt) 4 dE 



da 



(90) 
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which may be rewritten as 

^=i(^^J da ^ iaE J^ x U^ x ' a ' m) ' (91) 

with the "single particle integral" J; given as 

I(a,X,m) = [ d 3 p e iX P- ia ^ m2+p2 . (92) 



As proven in [25,26] and shown in Appendix A, one may 
write 



I(a, A, m) = 2-7T m 



(2) 

with E 2 representing a Hankel function. The phase 
space integral may thus be written as 



2 m 2 a itV) 



x 2 - A 2 2 



E^\m\J a 2 — A 2 ), (93) 



<f> = (f)(E, mi, ... , m„) 



2 n + l 7r 2n 

(2tt) 3 



u 

8 = 1 



rfa a n e iai5 



/•oo l, 2 n 

x y '' A (a2 _ Ay n^ ) ( m -v^) (94) 

We are now considering the "zero-mass limit", i.e. we 
calculate 

<f> n (E) := lim <f>(E, mi, ... , m n ). (95) 

In this case we may expand the Hankel function about 
the origin and keep only the first term. We obtain 



E^\mi\/a 2 - A 2 ) 



4i 



7rm 2 (a 2 — A 2 ) ' 
which may be substituted into eq. (94), to obtain 

23n-3j-n pco-ie 



(96) 



ME) 



da a n e lab E 2n , (97) 



with 



E, 



dX 



A 2 



(a 2 - A 2 ) 2 



(98) 



As shown in Appendix B, the integration can be done, 
and one obtains 

iw (4n-4)! 1 

#2n = " 2 4n-3 (2n - 1) ! (2n - 2) ! a 4 "" 3 ' ( "' ) 

The phase space integral is thus been given as 



ME) 



(4n-4)! 



7r 2-n 2 n (2n- l)!(2n-2)! 



fca-ie e' aE 

xl da- 



v 3n — 3 ' 



(100) 



(101) 



By choosing the contour in the upper half-plane, we ob- 
tain for the integral 



2« 



(iE) 



3n-4 



(3n-4)!' 



(102) 



and so the phase space integral, in the zero-mass limit, 
is given as 



M E ) = lim <f>(E, mi, ... , m n ) 

m, — ^0 



(4n-4)! 



E 3n ~ 4 (103) 



2«-i (2n- l)!(2n-2)!(3n-4)! 



This expression, eq. (103), can be evaluated easily, and 
is the basis for calculating multiplicity distributions, as 
discussed in the next section. 



IX. MULTIPLICITY SPECTRA IN THE 
ZERO-MASS LIMIT 

In this section we demonstrate how to calculate mul- 
tiplicity spectra in the zero-mass limit, and compare the 
results with the outcome of our Monte Carlo procedure, 
introduced earlier. This is not only a valuable check of 
the complicated numerical procedures, but also a very 
useful tool for optimizing our algorithm [27]. 

In our statistical treatment, the weight for a hadron 
configuration K = {hi, . . . , h n } is proportional to the 
partition function Sl(A'); correspondingly, the probability 
P n to find n hadrons, is given as 



Pn 



- E 



(104) 



Here, s = |<S| is the number of hadrons in the basic had- 
ron set S = {a i, . . . , a s }, and K ni „ s is the configuration 
with n v hadrons of species <r v . The condition Yjn v q v = Q 
accounts for flavour conservation, q v is the flavour vec- 
tor of hadron species v, and Q is the flavour vector of the 
droplet. In the zero-mass limit, eq. (104) may be written 
as 



1 

~Z 



Cvol <f>n 



E 



n 



9, v 



(105) 



where eqs. (2,3) have been used. The g v in eq. (105) 
have a different meaning than the in eq. (3): g v is the 
degeneracy of hadron species <r v . Z is a normalization 
factor. The prefactor C vo i and the phase space integral 
<f> n (in the zero-mass limit, see eq. (103)), do not depend 
on ni,..., n s , but only on n, and therefore appear in 
front of the summation symbol. We define 



"7 Cvol —7 1 

Z nl 



(106) 



12 



SUBATECH-95-05 



MARCH 1995 



HD-TVP-94-24 




FIG. 5. Medium size droplet: The multiplicity n versus 
the number of iterations for the hadron sets 5i (a), 53 (b), 
and 55 (c). The MC results, averaged over 200 iterations 
(solid lines) and 20 iterations (dotted), are compared with 
the analytical results for the average multiplicity (dashed). 



FIG. 6. Small size droplet: The multiplicity n versus the 
number of iterations for the hadron sets 5i (a), 53 (b), and 
55 (c). The MC results, averaged over 200 iterations (solid 
lines) and 20 iterations (dotted), are compared with the ana- 
lytical results for the average multiplicity (dashed). 



which is equal to P n for the case of = 1 and only one 
hadron species. In general, we have 

Pn = P° £ nlfl^. (107) 

n 1 ...n s V~\ 
En IJ q IJ =Q 



This is the final result for the multiplicity distribution in 
the zero-mass limit, which can be evaluated numerically 
as long as s is not too large, for example for our "test 
sets" Si, S3, or S5 . 

In figs. 5 and 6 we compare Monte Carlo (MC) results 
for the multiplicity with the "analytical results" for the 
average multiplicity, 
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FIG. 7. Medium size droplet: Multiplicity distributions 
for the hadron sets <Si (a), 53 (b), and Ss (c). MC results 
(solid lines) are compared with analytical results (dashed). 



FIG. 8. Small size droplet: Multiplicity distributions for 
the hadron sets Si (a), 53 (b), and Ss (c). MC results (solid 
lines) are compared with analytical results (dashed). 



<n>=J2 nP n, (108) 

n = l 

with P n from eq. (107). We consider a medium size 
droplet (E = 10 GeV and V = 10 fm 3 ) in fig. 5 and a 
small size droplet (E = 2 GeV and V = 2 fm 3 ) in fig. 
6; in both cases we have zero net flavour (Q = 0). We 
observe, indeed, that the MC results converge towards 
the analytical value. 



We now turn to multiplicity distributions. In figs. 7 
and 8 we compare Monte Carlo (MC) results, again for a 
medium size droplet (E = 10 GeV and V = 10 fm 3 ) and 
for a small size droplet (E = 2 GeV and V = 2 fm 3 ), with 
the corresponding "analytical results" , obtained from eq. 
(107). Also the average values Nmc and N ana are shown. 
The MC results are obtained from a single run per spec- 
trum (20000 iterations for the 10 GeV droplet and 200000 
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FIG. 9. Medium size droplet: Multiplicity distributions 
for the hadron sets 53 (a) and Ss (b). Exact analytical res- 
ults (solid lines) are compared with approximate treatments 
(dotted). 



FIG. 10. Small size droplet: Multiplicity distributions for 
the hadron sets 53 (a) and Ss (b). Exact analytical res- 
ults (solid lines) are compared with approximate treatments 
(dotted). 



for the 2 GeV droplet), which provides an accuracy of 
about 1 % for the 10 GeV case and of few % for the 2 
GeV case for the average multiplicities. 

For the larger set S7, and in particular for the realistic 
set 1S9, the exact expression cannot be handled. In this 
case we use an approximation, by neglecting flavour con- 
servation, i.e., we ignore the condition Hn v q v = Q. In 
this case, using the obvious identity 



\v = l 



E 



s 1 s 

n^n 



we get 



, E 



(109) 



(110) 
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FIG. 11. Medium size droplet: Multiplicity distributions 
for the hadron sets 57 (a) and Sg (b). MC results (solid 
lines) are compared with the approximate analytical ones 
(dashed). 



which can be evaluated also for 1S7 and S9. But first we 
compare in figs. 9 and 10 the exact and approximate res- 
ults for the sets 1S3 and S5 . Whereas for the medium size 
droplet the difference is quite small, we observe some dis- 
agreement for the small size droplet. We now turn to the 
large hadron sets: In figs. 11 and 12, we plot multiplicity 
distributions for the sets 1S7 and S9, comparing MC res- 
ults with the approximate analytical spectra. The MC 



FIG. 12. Small size droplet: Multiplicity distributions for 
the hadron sets £7 (a) and £9 (b). MC results (solid lines) 
are compared with the approximate analytical ones (dashed). 



spectra are shifted towards somewhat smaller multipli- 
cities, which is consistent with the observation in figs. 9 
and 10 that the exact results are "left-shifted" compared 
to the approximate ones. 

The multiplicity distribution P n refers to "total mul- 
tiplicities" , counting all hadrons. More information 
provide so-called "partial multiplicities" , where only the 
hadrons of a certain species are counted. We therefore 
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FIG. 13. Medium size droplet: Multiplicity distributions 
for 7r + (a), 7r - (b), p (c), and p (d) for the hadron set S5. 
MC results (solid lines) are compared with the analytical 
ones (dashed). 



FIG. 14. Small size droplet: Multiplicity distributions for 
7r + (a), 7r - (b), p (c), and p (d) for the hadron set S5. MC 
results (solid lines) are compared with the analytical ones 
(dashed). 



introduce the multiplicity distribution of hadron species Q ^-^ -JL gj 1 " 



jj, as 



=E p ° e « ! n vr- ( ii2 ) 



p k = \ E ^»...».). ( m ) 



a l • • - 1 n ii + l • • - n s V~\ 



"1 v iV+i'"*' j n order to achieve formal similarities to earlier formulas, 



where is fixed, and all other multiplicities n v (with 
v 7^ \i) are summed over. We may write 



we rewrite eq. (112) as 
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FIG. 15. Medium size droplet: Multiplicity distributions 
for 7r + (a), (b), p (c), and A (d) for the hadron set Sg. 
MC results (solid lines) are compared with the approximate 
analytical ones (dashed). 



FIG. 16. Small size droplet: Multiplicity distributions for 
7r + (a), (b), p (c), and A (d) for the hadron set Sg. 

MC results (solid lines) are compared with the approximate 
analytical ones (dashed). 



E 



9^ P° 



(113) 



e (»-».)! n 

u = l 



n 11 

M ilL. 



where and implies looping over v, except v = 

[i. The expression {} in eq. (113) has the same structure 
as the summation in eq. (107) for total multiplicities, 
so the same numerical procedures may be employed. In 
figs. 13 and 14, we show some some multiplicity spectra 
for specific hadrons. We compare Monte Carlo (MC) 
results, again for a medium size droplet (E = 10 GeV 
and V = 10 fm 3 ) and for a small size droplet (E = 2 
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GeV and V = 2 fm 3 ), with the corresponding "analytical 
results", obtained from eq. (113). 

For the larger sets S7 and S9, we use the approxima- 
tion of neclecting flavour conservation. In this case, the 
simplification eq. (109) may be used, and we obtain 



E 



9?* Pn 



(114) 



In figs. 15 and 16, we show some multiplicity spectra for 
specific hadrons for the large hadron set S9. We com- 
pare Monte Carlo (MC) results, again for a medium size 
droplet (E = 10 GeV and V = 10 fm 3 ) and for a small 
size droplet (E = 2 GeV and V = 2 fm 3 ), with the corres- 
ponding "approximate analytical results" , obtained from 
eq. (114). 



APPENDIX A 

In the following, we prove the relation 
I = J d 3 p e 8 '^P- 8 '«V m2 +P 2 

2_2 a rr(2) 



x 2 - A 2 2 



H^\m\/ a' 2 — A 2 ), (Al) 



(2) 

with H 2 representing a Hankel function, following 
[25,26]. 

Introducing polar coordinates , we have 
I = 2ir [ dpp 2 [ dcos^e iXpcos,} e^^+P 2 (A2) 



Jo J 

f°° 1 1 

= 27r/ dpp 2 —(e iXp -e- iXp )e- ia V m2+p2 (A3) 

n /■CO 

= - / dppe iXp - ia v m2 +P 2 . (A4) 
l * J-00 

We define £, A m , and a m via 

p=msinhC, A m = Am, a m = am, (A5) 
to obtain 

I=— rfCm 2 cosC sinCe iA ™ smhC - ia "« coshc (A6) 
*A ./_„ 



27rm 3 d r °° 
A m dX m 7_oq 



d( cos C e 



z'A m sinh £ — z'a m cosh £ 



We introduce an auxiliary variable <p via 



sinh ip 



cosh ijC 



/n- 2 _ \2 ' r Ay 2 _ \2 

V m m V m m 

The exponent in eq. (A7) may thus be written as 



(A7) 



(A8) 



i\ m sinh ( — ia m cosh ( = 

= "«V< - A2, (A9) 
x (— sinh ijC sinh C + cosh y> cosh £) 

= -i v /al-\l cosh(C-^). (A10) 



(All) 



With the integration variable 

using the identity 

cosh C = cosh ^ cosh <,£> + sinh £ sinh (A12) 
eq. (A7) may be written as 
j 2irm 3 d a m 



A m d\ m ^/a^-X^ 
x / ^cosh^e-V«™+ A ™ cosh e. (A13) 



Introducing 



A m 



m ' r/X r r/r ' 



(A14) 



7 = 27rm 3 a m — / d£ cosh £ e"""* coshe . 



z m dz m z m j _ qq 



(A15) 

Using the identities (see [25,26]) 

rf^COsh^ e-^ C ° She = -TTfrp^Ze''*) = -TTffp^z) 

) 

(A16) 



z dz z z l 



(A17) 



with ifj 8 - 1 representing Hankel functions, we obtain 

1 v( 2 ) 



7 = 27r 2 m 3 a m — H^\z m ). 



(A18) 



Using the definitions for z m , a m , and A m , we obtain the 
desired identity, eq. (Al). 



APPENDIX B 



We are going to prove the relation 
A 2 



H, 



dX- 



(a 2 - A 2 ) 2 " 
iw (4n-4)! 



2 4n-3 (2n - 1)! (2n - 2)! a 4 "" 3 ' 



(Bl) 
(B2) 
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for a = Rea — is. The integrand may be written as 

A 2 1 1 (m] 

(a 2 -A 2 ) 2 " a (a 2 -A 2 ) 2 " (a 2 - A 2 ) 2 "" 1 1 ' 

and, correspondingly, i?2n may be expressed as 

H 2n = a 2 A 2n - A 2n -i, (B4) 
with A n being defined as 

1 



(a 2 -\ 2 y 



The integrand has one pole in the upper half-plane, at 
A = —a = —Rea + is, since we are considering the case 
a = Rea — is. So we have 

A m = 2«Res {(a 2 - A 2 )" m } A= _ a . (B6) 
Expanding (a 2 — A 2 ) _m around s = A + a, we obtain 



(a 2 - A )~ m = (2a)" m (1 )" m £ 

2a 



s 

"2^ 



= (2a)" m £- m ^ 

i 

We can read off the residue of (a 2 — A 2 ) _m and obtain 



An 



27ri(2a)" 



(-2a) 1 



y m — 1, 
2«(2a) 1 - 2m (-l) 1 - m 

— m(— m — 1) . . . (—2m + 2) 



2wi (2a) 



12 . . .(m - 1) 
i- 2m (2m -2)! 



(m — 1)! (m — 1)! 
Using eq. (87), we obtain 

,(2a) 1 " 4n (4n-2)! 



#3 



2wi 



(2n- l)!(2n- 1)! 
(2a) 3 " 4n (4n-4)!" 
'(2n-2)!(2n-2)!_ 
— 7tz (4n-4)! 
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